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ABSTRACT 

We combine the physics of the ellipsoidal collapse model with the excursion set the- 
ory to study the shapes of dark matter halos. In particular, we develop an analytic 
approximation to the nonlinear evolution that is more accurate than the Zeldovich 
approximation; we introduce a planar representation of halo axis ratios, which allows 
a concise and intuitive description of the dynamics of collapsing regions and allows 
one to relate the final shape of a halo to its initial shape; we provide simple physi- 
cal explanations for some empirical fitting formulae obtained from numerical studies. 
Comparison with simulations is challenging, as there is no agreement about how to de- 
fine a non-spherical gravitationally bound object. Nevertheless, we find that our model 
matches the conditional minor-to-intermediate axis ratio distribution rather well, al- 
though it disagrees with the numerical results in reproducing the minor-to-major axis 
ratio distribution. In particular, the mass dependence of the minor-to-major axis dis- 
tribution appears to be the opposite to what is found in many previous numerical 
studies, where low-mass halos are preferentially more spherical than high-mass halos. 
In our model, the high-mass halos are predicted to be more spherical, consistent with 
results based on a more recent and elaborate halo finding algorithm, and with obser- 
vations of the mass dependence of the shapes of early- type galaxies. We suggest that 
some of the disagreement with some previous numerical studies may be alleviated if 
we consider only isolated halos. 

Key words: galaxies: clustering — cosmology: theory — ellipsoidal collapse, dark 
matter. 



1 INTRODUCTION 

A snapshot of a high-resolution dark matter simulation, 
at relatively low redshift, reveals that halos are neither 
spherically symmetric nor smooth. This is because, with- 
out even considering all the small scale complications arising 
from baryonic physics (i.e. pressure effects, merging, cooling, 
heating), the statistics of a Gaussian random field implies 
that spherically symmetric initial configurations should be a 
set of measure zero (Doroshkevich 1970). While it has long 
been recognized that spherical collapse is too idealized to be 
realistic, and that the true gravitational process must, at the 
very least, be ellipsoidal (Icke 1973; White & Silk 1979; Bar- 
row & Silk 1981; Kuhlman et al. 1996), most analytic mod- 
els of structure formation make the simplifying assumption 
that gravitationally bound objects are spherical, and formed 
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from a spherical collapse (i.e. Gunn & Gott 1972; Lahav et 
al. 1991; Lacey & Cole 1993; Lokas & Hoffman 2001). 

It has also become common practice to identify halos in 
simulations using a spherical overdensity algorithm, which 
finds the mass around isolated peaks in the density field such 
that the mean interior density is A times the background 
density, where the value of A is motivated by the spherical 
top-hat model (Lacey & Cole 1994). This is despite the fact 
that dark matter halos that form in simulations are rather 
elongated, and, in general, strongly triaxial: close to prolate 
(minor and intermediate axes are comparable in size to each 
other, and much smaller than the biggest axis) in the central 
parts, and rounder in the outskirts (Barnes & Efstathiou 
1987; Frenk et al. 1988; Dubinski & Carlberg 1991; Katz 
1991; Warren et al. 1992; Jing et al. 1995; Thomas et al. 
1998; Jing & Suto 2002; Allgood et al. 2006; Bett et al. 
2007; Diemand, Kuhlen & Madau 2007; Hayashi et al. 2007; 
Kuhlen et al. 2007; Munoz-Cuartas et al. 2010; Wang et al. 
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2010). This shape variation leads to a significant variance in 
local properties, compared to the spherically averaged value 
at a given radius, that are now becoming of interest. 

Studying and quantifying the degree of halo triaxiality 
is of broader importance. In fact, in the current paradigm 
of hierarchical clustering, dark matter halos are the hosts 
within which gas cools and collapses to form galaxies (White 
& Rees 1978; White & Frenk 1991), thus making them the 
building blocks of the large scale structure (LSS) of the 
Universe (Cooray & Sheth 2002). Hence, understanding the 
assembly histories, kinematics, clustering, and fundamental 
structural properties of halos - such as their intrinsic shapes 
- is the first necessary step in understanding the properties 
of galaxies (Mo, Mao & White 1998; Dutton et al. 2007; 
Diemand et al. 2008). In turn, the formation of dark mat- 
ter halos affects the properties of the galaxies hosted by the 
halos; therefore, inspection of the galaxy distribution in red- 
shift surveys such as the SDSS (York et al. 2000) allows one 
to relate the properties of galaxies to those of their host ha- 
los. Moreover, the characteristic density of a halo appears 
to track the mean density of the Universe at the time of 
its formation (e.g. Zhao et al. 2009), leading to a quasi- 
universal profile (Navarro, Frenk & White 1996; Kormendy 
& Freeman 2004), although self-similarity is not preserved: 
different halos cannot be rescaled to look alike (Navarro et 
al. 2004, 2010; Merritt et al. 2006; Gao et al. 2008; Reed et 
al. 2010). 

Triaxiality also has a number of observationally rele- 
vant implications. For example, modeling of dark matter 
halos beyond the spherical approximation is crucial in un- 
derstanding the nonlinear clustering of halos and dark mat- 
ter (Sheth, Mo & Tormen 2001), the formation and evolution 
of galaxies (Cole & Lacey 1996), and their relation to the 
cosmic web (Shen et al. 2006; Hahn et al. 2007; Lee, Hahn 
& Porciani 2009a,b; Forero- Romero et al. 2009; Pogosyan et 
al. 2009; Sousbie et al. 2009; Park, Kim & Park 2010). In 
particular, the higher order statistics of the nonlinear den- 
sity field is sensitive to halo triaxiality (Smith & Watts 2005; 
Smith, Watts & Sheth 2006). 

Triaxiality represents a useful framework for the non- 
spherical modelling of the intracluster gas, which recent ob- 
servations suggest will be key in deriving more accurate tem- 
perature profiles of X-ray clusters, and in general for cosmo- 
logical parameter determinations via the Sunyaev-Zeldovich 
effect (Lee et al. 2005). For example, Kawahara (2010) has 
derived the axis ratio distribution of X-ray clusters in the 
XMM-Newton catalog of Snowden et al. (2008) and con- 
firmed that the typical X-ray halo is well approximated 
by a triaxial ellipsoid. And recently, Morandi, Pedersen & 
Limousin (2010) have presented the first determination of 
the intrinsic triaxial shapes and three-dimensional physi- 
cal parameters of both dark matter and the intra-cluster 
medium for the galaxy cluster Abell 1689. 

Triaxiality is also useful in predicting - and hence can be 
constrained by - a variety of gravitational lensing observa- 
tions, including weak and strong lens statistics (Bartelmann 
1995; Van Waerbeke et al. 2000; Schulz et al. 2005; Bradac 
et al. 2006; Carbone et al. 2006; Bernstein 2007; Riquelme & 
Spergel 2007; Broadhurst et al. 2008; Limousin et al. 2008; 
Schneider & Er 2008; Mandelbaum et al. 2008, 2009; Zitrin 
et al. 2009; Bernstein & Nakajima 2009; Jonsson et al. 2010), 
and gravitational flexion (e.g., Hawken & Bridle 2009). By 



measuring the shapes of dark matter halos, galaxy-galaxy 
lensing can provide constraints on galaxy formation models 
and the nature of dark matter (Hoekstra et al. 2004; Man- 
delbaum et al. 2006; Parker et al. 2007). 

Surveys like ESA's Euclid mission 
(http://sci.esa.int/euclid) will in fact provide accurate 
data for shape estimates through "cosmic shear", a direct 
measure of the metric fluctuations in the Universe (Hoekstra 
& Jain 2008; Bernstein 2010; Rhodes et al. 2010), which in 
turn constrain dark energy properties (Albrecht et al. 2006). 
A primary source of noise in such measurements is due 
to the difficulty in distinguishing between intrinsic galaxy 
shapes and shape distortion due to lensing (Bartelmann 
& Schneider 2001; Refregier 2003; Hoekstra et al. 2005; 
Mandelbaum et al 2006; Bridle et al. 2009). Hence accurate 
modelling of the correlated shapes and orientations dark 
matter halos can be extremely useful. The higher order 
statistics of the nonlinear density field in such surveys is 
also sensitive to halo triaxiality (Smith et al. 2006). 

On galaxy mass scales, an understanding of halo tri- 
axiality provides useful input to studies of galactic disks in 
triaxial halos. E.g., Jeon, Kim & Ann (2009) considered the 
fundamental dynamics between the disk and the axisym- 
metric or triaxial halo, and Valluri et al. (2010) analyzed 
the orbital structure of dark matter particles in TV-body 
simulations in an effort to understand what is the phys- 
ical mechanism driving shape changes caused by growing 
central masses (also see Debattista et al. 2008). This shape 
change reconciles the strongly prolate-triaxial shapes found 
in collisionless iV-body simulations with observations, which 
generally find much rounder halos (see for example Banerjee 
& Jog 2008). Finally, resolving the fine grained structure of 
galaxy mass halos enables one to make more realistic pre- 
dictions for direct and indirect dark matter detection exper- 
iments (e.g. Giocoli, Pieri & Tormen 2008). 

The shapes of dark matter halos can be quantified using 
high-resolution JV-body simulations of hierarchical gravita- 
tional clustering. This is currently the best way to address 
many of the tasks above. Simulations are becoming increas- 
ingly accurate in mass and spatial resolution: e.g., the Mil- 
lennium Run (Springel et al. 2005), the Horizon Run (Kim 
et al. 2009), the Millennium II (Boylan-Kolchin et al. 2009) 
and the Bolshoi Simulation (Klypin et al. 2010). 

In these numerical studies, many different aspects have 
been considered over a wide range of physical scales and cos- 
mic histories. Common findings suggest that the total mass 
is the key in determining the final shapes of halos, although 
other environmental parameters may play a role in the pro- 
cess. In general, halos do exhibit a rich variety of shapes 
with a preference for prolateness over oblateness. More mas- 
sive halos tend to be less spherical and more prolate, and 
are preferentially aligned with primordial filaments, while 
less massive halos are in general rounder (e.g., Jing & Suto 
2002; Allgood et al. 2006; Araya-Melo et al. 2009). On the 
other hand, perturbation theory suggests that the most mas- 
sive objects should be spherical (Bernardeau 1994; Pogosyan 
et al. 1998), and Lemson (1995) showed that the spherical 
model does provide a good description of the evolution of 
the spherically averaged profile. More recently, Dalai et al. 
(2008) reported that massive halos are indeed very well de- 
scribed by the spherical collapse model. This was recently 
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confirmed by Park et al. (2010), who provided a clear expla- 
nation for the discrepancy with previous measurements. 

The shapes of subhalos are similar to those of host ha- 
los, but subhalos tend to be a bit rounder, especially the 
ones near the host halo center. Tidal interactions make in- 
dividual subhalos rounder over time and they tend to align 
their major axis towards the center of the host halo. Forma- 
tion of halos is also affected by the large-scale environment, 
which may have an impact on their shapes, and those shapes 
can be modified by galaxy formation as well. Subhalos are 
not the subject of our study. 

There are many more numerical studies of halo shapes 
than analytic models. This is because the formation, evo- 
lution and virialization of dark matter halos is complex; no 
rigorous analytic techniques are available for use in both the 
linear and the nonlinear regimes. In addition, choosing the 
appropriate definition of halo shapes is subtle. For exam- 
ple, Eisenstein & Loeb (1995) describe an analysis of halo 
shapes which uses the ellipsoidal collapse model of Bond & 
Myers (1996). However, as we discuss below, their definition 
of a halo differs from the more commonly accepted defi- 
nition. Moreover, they present results for collapsed objects 
that had the same initial density. Since the time it takes for 
to collapse is a complicated function of density and shape 
(Bond & Myers 1996), this means that they compare objects 
of one shape at one time with those of another shape at a 
different time. This is rarely measured in simulations: the 
shape distribution of most interest is, of course, that for a 
fixed time (e.g., halos at z = 0). 

The main purpose of the present work is to provide 
a simple model for the distribution of halo shapes at any 
given time as seen in the simulations, starting from first 
principles. Following Rossi (2008), our analytic prescription 
has two independent parts: the first is a scheme for how an 
initially spherical patch evolves and virializes; the second is 
the correct assignment of initial shapes to halos of different 
masses. 

In this respect, our model is similar in philosophy to 
that of Lee, Jing & Suto (2005). However, there are impor- 
tant differences. (1) They assume the Zeldovich approxima- 
tion remains valid even during the nonlinear regime, where 
it is known to fail; (2) They assume a spherical collapse 
threshold for the formation of halos, or an empirical recipe 
based on Lee & Shandarin (1998). In contrast, because we 
are modelling triaxial objects, we self-consistently use ellip- 
soidal, rather than spherical collapse dynamics to generate 
our predictions. 

To describe the evolution of non-spherical structures we 
adopt the ellipsoidal collapse model of Bond & Myers (1996), 
which was used by Sheth, Mo & Tormen (2001) to estimate 
of how the abundance of dark matter halos depends on halo 
mass. In this model, dark matter halos are identified with 
ellipsoids which have collapsed completely along all three 
axes (we show below that, in effect, the Eisenstein & Loeb 
1995 definition corresponds to collapse along just two axes). 
In this framework, the time required to collapse depends on 
the overdensity Si and size Ri of the initial patch, and on 
the surrounding shear field, parametrized by its ellipticity 
e and prolateness p. Requiring the collapse to happen at a 
given time makes Si a function of e and p (Sheth et al. 2001): 
S cc (e,p). The combination of Si = S cc , e and p determines 
the axis ratios of the object at all times, and, in particular, 



at the final time. Thus (e,p) establishes the time of collapse 
(it was this fact which was exploited by Sheth et al. 2001), 
as well as the axis ratios at collapse (a fact we exploit here). 

The second part is the correct assignment of (e,p) values 
to halos of different masses. We do this following Sheth & 
Tormen (2002) (also see Chiueh & Lee 2001; Sandvik et al. 
2007). In essence, the correct (e,p) distribution is specified 
by the statistics of Gaussian random fields. In a Gaussian 
random field, the distribution of (e,p) values depends on the 
size and overdensity of the patch: g(e,p\S, R). Massive halos 
form from larger patches in the initial conditions than do less 
massive halos, so we expect the distribution of initial (e,p) 
values, and hence the distribution of final axis ratios, to also 
depend on halo mass. Thus, although the generic evolution is 
initially towards an oblate, pancake-like structure, followed 
by a shift towards a more prolate shape, as the other axes 
also begin to shrink, quantifying this evolution, and merging 
it with the correct (mass dependent) initial distribution of 
shapes is the main focus of this work. 

The outline of this paper is as follows. In Section [2] we 
present the ellipsoidal model for the gravitational collapse of 
an initially spherical patch; we discuss a reasonably accurate 
analytic approximation to the evolution, more rigorous than 
the Zeldovich approximation; we introduce the axis ratio 
plane, which provides a concise description of the dynamics 
of collapsing regions and allows one to relate the final shape 
of a halo with its original, pre-collapsed, shape. In Section 
Owe present the full model for halo shapes, we explain how 
the initial conditions are obtained via the excursion set al- 
gorithm, and we expand on the prolateness distribution - 
crucial in understanding our main results. In Section U the 
model is contrasted with high resolution iV-body simula- 
tions. A simple explanation of an empirical relation found 
by Jing & Suto (2002) is given, and a number of caveats in 
the comparison model/simulations are highlighted. Finally, 
Section [S] discusses in detail limitations in the modelling and 
difficulties arising from numerical studies, and suggests fu- 
ture improvements. 

In our calculations we assume a spatially flat cosmolog- 
ical model with (Qm,Qa, h) = (0.3,0.7,0.7), where Qm and 
Qa are the present day densities of matter and cosmological 
constant scaled to the critical density. We write the Hub- 
ble constant as Ho = lOO/i km s" 1 Mpc -1 . Regarding the 
overall notation, we always use the subscript i to denote ini- 
tial quantities, and the subscript / for final quantities. The 
expansion factor of the Universe, with mean density p, is 
represented by the small case letter a - not to be confused 
with the capital letter A, which denotes the axis ratios of an 
ellipsoidal patch. 



2 NONLINEAR DYNAMICS 

The nonlinear gravitational evolution of a medium is com- 
plex, even without complications arising from gas dynamics. 
It involves smooth accretion, tidal interactions with the en- 
vironment as well as violent episodes of collisions with other 
halos, merging and fragmentation. In what follows, we first 
briefly summarize the key aspects of the ellipsoidal collapse. 
We then discuss an analytic approximation for the evolu- 
tion, and introduce a planar representation of halo axis ra- 
tios; this allows for a concise description of the dynamics 
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of collapsing regions, and provides a mapping between the 
initial and final shape of a halo. 



2.1 A model for the gravitational collapse 

The simplest description of the gravitational evolution and 
virialization of a cosmic structure is the spherical collapse 
model (Gunn & Gott 1972; Gott 1975; Gunn 1977; Fill- 
more & Goldreich 1984; Bertschinger 1985; Mo & White 
1996). In this framework, an isolated overdensity in an oth- 
erwise unperturbed universe first expands with the Hubble 
flow, then turns around and collapses. However, the initial 
shear field, rather than the density, has been shown to play 
a crucial role in the formation of nonlinear structures (Zel- 
dovich 1970; Hoffman 1986, 1988; Peebles 1990; Dubinski 
1992; van de Weygaert & Babul 1994; Audit & Alimi 1996; 
Audit, Teyssier & Alimi 1997). Therefore, refinements to the 
spherical approximation which include local shear effects can 
be obtained by introducing an ellipsoidal collapse scheme 
(Bond & Myers 1996; Eisensein & Loeb 1995; Monaco 1995, 
1997, 1998; Lee & Shandarin 1998; Chiueh & Lee 2001; 
Sheth, Mo & Tormen 2001; Sheth & Tormen 2002; Ohta 
et al. 2004; Shen et al. 2006; Sandvik et al. 2007; Desjacques 
2008; Desjacques & Smith 2008; Rossi 2008). 

In this study we adopt the homogeneous ellipsoidal 
model in the form of Bond & Myers (1996), although the 
ellipsoidal collapse has a long history (Lin, Mestel & Shu 
1965; Icke 1973; White & Silk 1979; Barrow & Silk 1981). To 
lowest order, their algorithm reduces to Zeldovich's (1970) 
approximation, and so linear theory is reproduced. In this 
framework, an initially spherical patch of initial size Ri with 
overdensity Si is distorted by the shear field into a collaps- 
ing homogeneous ellipsoid. The exterior tidal force arising 
from the matter outside of the ellipsoid is completely deter- 
mined by the volume-averaged strain of the ellipsoid. The 
details of the substructure in the interior are ignored, and 
the strongly nonlinear internal dynamics of the collapsing 
patch are assumed to be largely decoupled from the weakly 
nonlinear dynamics describing the motion of the patch it- 
self; in this respect, the homogeneous ellipsoid picture may 
be thought of as a tensor virial theorem approach to the 
average interior dynamics. In fact, in the nonlinear regime 
the one-to-one correspondence between the external tidal 
field and the local strain tensor is no longer true on small 
scales where the rms density fluctuation a 3> 1; therefore, 
one would not expect the simple ellipsoidal model to apply 
in this regime. 

There is another sense in which this approach is only a 
simple approximation. It assumes that the inertia tensor of 
the final bound object is perfectly correlated with the local 
tidal tensor in the initial Lagrangian space. Measurements 
in simulations show that the two tensors are not perfectly 
correlated (e.g. Lee & Pen 2000) . On the other hand, the cor- 
relation is stronger than naive tidal torque theories predict 
(e.g., Porciani, Dekel & Hoffman 2002), so the assumption 
of perfect correlation is a useful idealization. 

It is usual to characterize the initial shear field by the 
ellipticity e and prolateness p associated with the potential 
rather than with the density field (i.e. Bardeen et al. 1986). 
This is because the components of the 3x3 strain tensor 
are the second derivatives of the potential. The eigenvalues 
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Figure 1. Evolution of axis lengths in our model, in physical 
units. The times at which different axes freeze out are determined 
by the initial (e,p, 8) values, as specified in the panels. 



of the initial strain tensor are related to the initial density 
contrast and to the shear ellipticity and prolateness by: 

Ai(*i) = ^(l + 3e + p), (1) 
A 2 (*0 = ^(l-2p), (2) 

A3(ti) = ^p(l-3e + p). (3) 

Note that ^\ Aj = <5j. If S > 0, then e ^ and — e ^ p ^ e, 
so Ai > A2 > A3. 

If we denote with A k the scale factors for the three prin- 
cipal axes of the ellipsoid, then the initial conditions are set 
by the Zeldovich approximation, both for the displacement 
and the velocity fields: 



A k (ti) = o(ti)[l-A k (ii)] 
dA k (ti 



dt 



A k (ti) - o(ti)A k (ti) 



d \nD 



d lna 



(4) 
(5) 



Notice that A3 ^ An ^ Ai. The subsequent evolution is 
given by 



dt 2 



k = Q A H§A k - 4nGpA k 



1 + 5 b' k 5 , ' 

— ^ 1 2 ' A ext,k 



(6) 
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where p is the mean density of the Universe, 8 the relative 
overdensity, and b' k = bk — 2/3 and A[, xt k account for the 
interior and exterior tidal forces. In particular, 



fck(t) 



and 



[n 



A m (t) 



lAt(t)+T]ilUAf(t) + r]^ 



(7) 



A e xt,k(0 



D(t) 



D(U) 



A k (ti) - <5(t ; )/3 



(8) 



with D the linear theory growing mode. Other possibilities 
for Ag Xt k include the 'nonlinear' (Bond & Myers 1996), or 
the 'hybrid' (Angrick & Bartelmann 2010) approximations. 
However, in our case we are concerned with later times, 
where a different choice for the external shear field does not 
significantly affect our conclusions. 

If e = p — 0, then all three eigenvalues equal Si/3, hence 
all three axes have the same length initially. In this case, the 
exterior anisotropic tidal force is zero and b' k = 0, so one gets 
the usual cycloid solution for a closed universe. Hence, all 
three axes evolve similarly, so the object remains spherical, 
and the time to collapse is determined by one number: Si. 
But for more general initial values of 8(ti), e and p, and an 
initial redshift, equation (|6]) must be solved numerically for 
each axis A k . Generically, a triaxial object has three critical 
times, corresponding to the collapse along each of the three 
axes. In this case, the shortest axis, Ai, collapses first and 
A3 collapses last. 

Figure [T] shows how the lengths of the three axes evolve 
in the model, in physical units. In all cases we set 5(ti) = 
[D(ti)/D(to)]S(t ), with 5(t ) = 1.6753 being the usual crit- 
ical value associated with spherical collapse in the adopted 
cosmology at z = 0. Our numerical calculations start at a 
time ti, which corresponds to a redshift Zi = 39. Each panel 
show results for a different pair of (e,p). For a given e, p > 
implies a pancake-like structure (i.e. one short axis and two 
long), while p < results in filament-like structures. The 
main point of the figure is to illustrate that, in this model, 
a given pair (e,p) determines the axis ratios of the object at 
all times, and, in particular, at the final time. 

Halos are identified with objects that have collapsed 
along all three axes. Bond & Myers (1996) stop collapse 
along axis k by simply freezing A k once a critical radius 
^4cq,k = a/ r is reached during the infall phase, where a is 
the expansion factor of the Universe and typically the radial 
freeze-out factor f r = 0.177 is chosen in order to reproduce 
the 'viriaP density contrast of A = a 3 /(Ai A 2 A :j ) = 179 
familiar from spherical top-hat calculations in an Einstein- 
de Sitter model; f r must be computed for more complicated 
cosmologies. 

In our study, we slightly modify the stop criterion pro- 
posed by Bond & Myers (1996) as follows: we still freeze 
the value of axis i once a critical radius is reached during 
the infall phase (the radial freeze-out factor being chosen in 
order to reproduce the correct spherical virial density con- 
trast in the assumed cosmological model), but we progress 
the evolution in time till we reach the point in which each 
axis has collapsed completely, in order to consider fully re- 
laxed halos. Bond & Myers (1996) found that the time at 
which the longest axis of the ellipsoid freezes out is relatively 
insensitive to the exact value of / r , within a given cosmo- 
logical model. In this respect, a more sophisticated collapse 
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Figure 2. Comoving evolution of axis lengths in our ellipsoidal 
model. Solid lines are obtained by solving equation J6]l numeri- 
cally; dashed-dotted lines show equation H01 . 



criterion such as the one recently proposed by Angrik & 
Bartelmann (2010), and based on the tensor- virial theorem, 
does not affect our analysis. 

However, this stop criterion is rather different from that 
used by Eisenstein & Loeb (1995), and this difference does 
matter. In their prescription, the collapse is stopped at the 
time when a sphere, whose overdensity was that of the initial 
ellipsoid, would have shrunk to zero radius. This happens to 
be very close to the time when the intermediate axis col- 
lapses (Shen et al. 2006, and see discussion below), so it can 
be substantially before the time that the third axis collapses. 

Before moving on, we note that the critical density re- 
quired for collapse by the present time is well-approximated 
by 

Me,p)~ f ac (9) 

1-/V5(e 2 ±p 2 ) 

with P = 0.365. This will be useful in what follows. 



2.2 An analytic approximation 

To characterize the dynamics of a collapsing object, one 
must solve numerically the coupled partial differential equa- 
tions for the ellipsoidal collapse, as shown in Figure [T] How- 
ever, there is a useful analytic approximation to the exact 
solution which provides considerable insight. By extending 
results in White & Silk (1979), Shen et al. (2006) show that 



A k (t) a(ti) 



A k (ti) a(t) 



1 



Pit) 
D{ti) 
A h (ti 



Ak(ti) 



D(t) S(ti) a c (t) 



(10) 



D(ti) 3 a(t) 

where Ah(ti) = 3/^. J 4/ 1 (ti) and a c (t) is the expansion 
factor of a universe with initial density contrast <5(ii) = 




Figure 3. The 'axis ratio plane' for dark matter halos. Left, central and right panels show p = and p = ± e/2, respectively, when 
e = 0.15. The initial redshift is z\ = 39, for an initially spherical patch with overdensity 5{t\) = [D(t{) / D(to)]8(to), where <5(to) = 1.6753. 
See the main text for more details. 



y\ Aj (ti) . This approximation to the full ellipsoidal collapse 
model can be thought of as correcting the Zeldovich (1970) 
approximation by the factor by which it is wrong for a sphere 
(Lam & Sheth 2008). 

To first order in S(ti), 



Ai, 3 (t) a(tj) ^ a c (t) S(U) 

Ai, 3 (ti) a(t) ~ a(t) 3 

A 2 {i) a(tj) _ ae(t) , S(ti) 
A 2 {ti) a(t) ~ a(t) 



(P±3e) 



1 + 



D(i) 0e (t) 



D(ti) a(t) 



+ 



-2p 



1 + 



D(t) a e (i) 



D(ti) a(t) 



(11) 



(12) 



(13) 



Notice that when p = 0, then 

A 2 (t) a e (t) 
A 2 (U) ~ a(ti) ' 

in this case, the second axis evolves exactly as in the spher- 
ical model, since for the spherical case e = p — 0, Ak(ti) = 
A(*i) = 5(ti)/3 and so A k (t) = A(t) -> A(U)[a e (t)/a(U)}. 
In general, the second axis evolves very similarly to that 
for a spherical model with the same initial overdensity 5(ti). 
In this respect, a spherical collapse can be roughly seen as 
an 'imperfect' ellipsoidal model, where the virialization is 
identified with the collapse of the second axis. 

Figure[2]shows how well the approximation works in the 
LCDM model assumed in this study. Solid lines show nu- 
merical solutions of equation © in comoving units; dashed- 
dotted lines show the analytic approximation (equation [TUJ). 
As evident from the figure, the approximation is rather good 
until the collapse of the third axis, as expected. 

2.3 The axis ratio plane: dark halo evolution 

The evolution of triaxial objects can be conveniently de- 
scribed by a two-dimensional 'axis ratio plane', showing the 
shortest-to- longest (A1/A3) versus intermediate-to- longest 
(A2/A3) axis ratios. Figure [3] shows what our ellipsoidal 
collapse model predicts for a few different combinations 
of prolateness and ellipticity, as indicated in the panels, 
when S(ti) = [D(ti)/D(t )]6(t ) and % = 39. Note that 



S(to) = 1.6753 is again the usual spherical collapse linear 
value in the concordance cosmology, at the present time. 

The axis ratio plane provides useful insights into the 
evolution of dark matter halos. The trajectory of a collapsing 
object in the plane is as follows: down and to the left, as the 
shortest axis shrinks more rapidly than the other two (first 
line); then further to the left and slightly upwards, after the 
shortest axis has frozen out, while the second axis continues 
to shrink faster than the longest axis (second line); finally 
upwards and to the right, after the second axis has also 
frozen out, so only the longest is shrinking (third line), until 
the third axis also freezes out, at which point the object is 
defined as being virialized. 

The slope of the first line is steeper when p = e/2, 
shallower when p = —e/2. The evolution proceeds further 
down the plane as e increases, at fixed p, and is confined 
in the region between the bisector of the plane and the line 
characterized by p = e/2. Hence, when p = e/2 and e > 0, 
the object is more likely to form a pancake (i.e. lower right 
region of the plane), while when p = —e/2 a filament is more 
likely to occur (i.e. lower left region). Also, for small values 
of ellipticity, the final shape does not depart significantly 
from the spherical case (upper right region of the plane). 

The fact that p determines the slope of the initial mo- 
tion in the plane, whereas e controls how far the initial col- 
lapse progresses, can be understood as follows. Initially, i.e., 
when Ai <C 1, 



Ai 

A 3 

A2 
A 3 



1-Ai 



25(ii)e, 



I-A3 
1 — A3 A3 



(14) 
(15) 



Therefore, for a given 5(ti), the ratio A1/A3 is a function of e 
only: if e increases then A\ I A3 decreases, and at fixed A\ I A3 
(or, equivalently, for a given e) the prolateness p increases 
as A2IA3 increases. This effect is indicated in Figure [3] by 
the circular arrows. 
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Ai(t) A 2 (t) l-5(ti)(l + 3e + p)/3 



Sphere 



Sphere 



Mt) 



5(ii)(l-2p)/3 



{l-^^CSe + p) 


i , D(t) a c 
± D(ij) a 


(()" 
(f) _ 




j\ . 2<5(t,) a(t) 
| 1+ 3 a a (t)P 


, , £ft) a e (t) 
1 "^(77 


}■ 



(16) 



In particular, when p = 

4»(t) 



D(t) \ a{t) 



(l+3e){l+*(fa)e[l-(l+ 5 ^) 



\M1 (17) 



D(ti)J de(t) 

Motivated by this, we have parametrized the evolution 
in this plane as follows. If /(e) denotes the distance along 
the first line, then, for a given p, 

f(e) = ^/(1-A 2 /A a y + (1-A 1 /A s )2, (18) 
and 



A- A - k{p) A-r b - 

The constraint that Ay/ A3 = 1 when A2/A3 
that 6=1 — k(p), and therefore: 

A 2 



f(e) 



(19) 
1 implies 

(20) 



We found k{p) = 1.858 when p = 0, fc(p) = 3.803 when 
p — e/2, and k(p) = 1.253 when p = — e/2 for the collapse 
of the shortest axis. For the turnaround of the second line, 
parametrized in the same fashion, k(p) = 1.084 when p = 0, 
k(p) = 1.253 when p = e/2 and k{p) = 1.034 at p = -e/2. 
The slope of the final part of the evolution (third line) , when 
Ai and Ai are both frozen, is simply given by Ai(tf)/A2(ts). 
Hence, it approaches unity (i.e. parallel to the bisector of the 
plane) when A1/A3 ~ A2/A3, and it is always lesser than 
unity otherwise, since Ai(t{) ^ Ai(ti). 



3 INITIAL AND EVOLVED HALO SHAPES 

The discussion above shows how the initial values of S, e,p 
determine the future evolution of the object, including its 
shape. Our next task is to determine how the initial distri- 
bution of e and p values depends on halo mass. We use the 
algorithm of Sheth & Tormen (2002) to do this. Briefly, this 
requires the generation of the joint distribution of 5,e,p as 
a function of smoothing scale Ri, finding the largest scale at 
which S ^ S cc (e,p), and then associating the values (S cc , e,p) 
to the patch from which a halo of mass M = pAivRf /3 
formed. Since this algorithm has been used by other au- 
thors since (Sandvik et al. 2007), we do not provide details 
here, but refer the reader to Appendix B in Sheth & Tormen 
(2002), and to Bond et al. (1991), Bower (1991) and Lacey 
& Cole (1993) for further details. 



3.1 Mass-dependent shapes from scale-dependent 
initial conditions 

Figure [4] illustrates how our model for halo shapes works. 
The top left panel shows the distribution of initial axis ratios 
for small patches (M — 10 13 /i _1 Mq) which the ellipsoidal 
collapse model predicts are destined to become low mass 
halos at the present time. The bottom left panel shows the 
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Figure 4. Distribution of initial (top panels) and final (bot- 
tom panels) axis ratios in our model, for halos of mass M = 
10 13 ft- 1 M (left) and lO 14 h~ 1 M (right). Initial conditions are 
specified by combining the statistics of Gaussian random fields 
with the ellipsoidal collapse model, following Sheth & Tormen 
(2002). 



distribution of final axis ratios at virialization (see Section 
[2j|. The panels on the right show the analogous quantities 
for more massive halos (M = 10 14 ft _1 Mq). 

The plot shows that halos are, in general, predicted to 
be triaxial, with a slight tendency to be more prolate than 
oblate. In addition, massive halos are predicted to be more 
spherical than less massive halos, both initially and finally. 
This is primarily a consequence of the fact that the larger 
patches from which massive halos form have, on average, 
smaller values of (e,p). The trend is consistent with the find- 
ings of Bernardeau (1994), who argued that in perturbation 
theory larger halos are expected to be rounder. 



3.2 The conditional prolateness distribution 

The joint distribution of Ai ^ A2 ^ A3 was derived by 
Doroshkevich (1970) for Gaussian random fields (see Lam, 
Sheth & Desjacques 2009 for an extension to non-Gaussian 
fields that are of the 'local' type). The equivalent expression, 
in terms of e, p and S, is: 

g(e,p,S\a) = p(5\a) g(e,p\5, a) 
_ s 2 

e 1125 , 2 _2^<5^_j4 C 3 e 2 +J . 2 ) , . 



Doroshkevich 's formula is the product of two independent 
distributions, a Gaussian for S/a, and one which is a com- 
bination of the other five independent elements of the de- 
formation tensor of which the Aj are the eigenvalues (Sheth 
& Tormen 2002) . Integration of equation 1)211) over 5 yields 
the joint distribution of ellipticity and prolateness, g(e,p\a). 
A further integration over p gives the distribution of ellip- 
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Figure 5. Distribution of prolatencss at fixed halo mass and ellipticity, g c (p\e,a), in our model (dotted histograms). Theoretical curves 
(solid lines) are from equation J22I) with different values of e, as specified in each panel. The initial distribution g(p\e, a, 8) is also shown, 
with long-dashed lines (i.e. equation [23} , as well as its approximation (i.e. dotted lines). See the main text for more details. 



ticities as a function of halo mass (<r), /(e|cr). This integral 
can be computed analytically (c.f. equation 24 of Lam et al. 
2009). 

For a given ellipticity, we can predict the distribution of 
p's that will produce collapsed objects at any given redshift 
(say, z = 0). In detail, for 5 ec (e,p) as in equation @, we 
define €>(<5, e,p) = 1 if 6 8 ec (e,p) and Q(S,e,p) = other- 
wise. Then we impose the requirement that the perturbation 
had collapsed (at z = 0) by computing 



g c (p\e,cr) = 



9c(e,p\ff) 



Jl e d P 9c(e,p\a) 
J °° dSg(e,p\S,a)p(S\a) @[5 > S ec (e,p)] 
Jl e d P / °° dS 9(e,p\S, a) p(S\a) Q[S > S cc (e,p)] 



(22) 



Plots of equation (]22[) are shown in Figure[5]with solid lines. 
We can further quantify if the difference between g(p\e, S, a) 
and g c (p\e,a) is relevant (i.e. dependence of barrier and 
epoch) by computing: 



g(p\e,8,a) 



g(e,p\S,a) 



Jl e A P g(e,p\S,a) 



25^(l-p 2 /e 2 )e^ 2 / 2 » 1 -P 2 / e2 )/e 

. i) e (5y 2 /2) er f( x /5y 2 /2) + Wy ' 



(23) 



where y = eS/cr. Long-dashed lines in Figure [5] show this 
distribution: as evident from the panels, the difference be- 
tween (1221) and (|23[) is negligible. In fact, if one ignores the 
exponential prefactor in equation (|2ip . then it is straight- 
forward to show that 

g(p\e, 6, a) ~ g(p\e) = (3/4) (1 - p 2 /e 2 )/e. (24) 

Dotted curves in Figure [5] show that this is a very good 
approximation. 

3.3 Universal conditional axis ratio distributions 
at late times 

In equation Q22p . the mass and epoch dependencies actually 
cancel out, i.e. g c (p\e,a) = g c (p\e). This suggests that non- 
linear evolution will not change the conditional distribution 



g c (p\e,a) from that which is set by the initial conditions. 
To see why, recall that the axis ratio A1/A3 is specified once 
S(ti) and the ellipticity are fixed (c.f. equation ll4|l . Moreover, 
if 6{ti) and e are fixed, then A1/A2 = (Ai/A 3 )/(A 2 /A 3 ) is 
only a function of the prolateness (c.f. equation 1 1 5 p . Hence, 
for a fixed mass, equation (I22C can be equally expressed in 
terms of (e,p), or as a function of (A1/A2, Ai/As). I.e., if 
we define A12 = A1/A2 and ^13 = A1/A3, then 



p(A 12 \Ai 3 ,5,a) = 



g(A 12 ,A 13 \S,a) 
p(A 13 \8, a) 
3 



2(1- 



1 



(2Ai 



Ax 



(1 



x exp{-^(2A 12 -l- J 4 13 ) 2 },(25) 
where the final expression follows from setting 



g(A 12 , A l:i \S,a) dAi2 dAi 3 = g(e,p\8, a) dei 



(26) 



At small masses, the exponential term — > 1, making this 
distribution almost universal, i.e. independent of mass and 
epoch. This remarkable feature will guide the interpretation 
of our results on halo shapes in the next section. 



4 COMPARISON WITH HIGH RESOLUTION 
SIMULATIONS 

The previous section showed how we combine the nonlinear 
ellipsoidal model for halo collapse and virialization (Section 
[2| with the excursion set theory (Section [21), to model the 
initial and evolved shape distributions of dark matter ha- 
los. Here, we investigate to what extent our simple analytic 
prescription is able to reproduce results from iV-body sim- 
ulations of Jing & Suto (2002). 

For their statistical analysis of halo shapes, Jing & Suto 
(2002) identified halos at z = 0, 0.5, 1 in simulations of a 
cosmology with (JIm = 0.3, JIa = 0.7, h = 0.7). The sim- 
ulations used -/V = 512 particles in a 100 h _1 Mpc box. 
Halos were identified using a friends-of-friends (FOF) algo- 
rithm with link- length b — O.ld, where d is the interparti- 
cle separation. Each of their FOF halos contains more than 
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Figure 6. Distribution of final axis ratios in our model for halos 
of different masses (solid contours), as specified in the panels. 
Each halo has been evolved using the ellipsoidal collapse model. 
The generic evolution is initially towards an oblate, pancake-like 
structure, followed by a shift towards a more prolate shape, as 
the other axes also begin to shrink. Dotted lines are the joint 
distribution which Jing & Suto (2002) report describe the halos in 
their simulations, in the concordance cosmology (Cm = 0.3, = 
0.7, h = 0.7). 



10 4 particles. The lower mass limit of their halo catalog is 
6.2 x 10 12 /i -1 Mq. They provided a set of fitting formulae 
which allow simple tests of our triaxial model: in particular 
the mass and redshift dependence of the axis ratios, and the 
probability distribution functions of the axis ratios. 

4.1 Joint distribution of axis ratios 

By averaging over the correct, mass dependent, distribu- 
tion of initial shape parameters, each evolved through the 
ellipsoidal collapse model, we are also able to make predic- 
tions for how the distribution of final shapes depends on halo 
mass. The solid contours in Figure [6] show the joint distri- 
bution of (A2/A3, A1/A3) we predict; contours show levels 
where the probability has fallen to 1/2, 1/4, 1/8, 1/16 of the 
maximum value. The dotted lines represent the joint distri- 
bution which Jing & Suto (2002) report describe the halos 
in their simulations. Note, for the moment, that our model 
appears to be in reasonable agreement with the simulations 
only around M ~ 10 13 h' 1 M e ~ M» . While this agreement 
is gratifying, the fact that the mismatch is mass-dependent 
is potentially troubling. Jing & Suto find that the low mass 
halos are actually slightly more spherical than massive halos. 
Although the effect is weak, it has since been confirmed by a 
long list of authors (see Section[l|. In contrast, in our model, 
it is the massive halos which are expected to be more spher- 
ical (Figure [4}. The trend we find is consistent with pertur- 
bation theory (Bernardeau 1994; Pogosyan et al. 1998; see 
Park et al. 2007, Dalai et al. 2008 and Park, Kim & Park 
2010 for more recent discussion of this point). 



To show this more clearly, Figure [7] compares the final 
distribution of A1/A3 in our model (dotted) with the distri- 
bution in Jing & Suto's simulations (solid) at z = 0. Again, 
note that the agreement is reasonable around M*, but that 
the predicted distribution is broader than the simulations at 
low and high masses. Note also that the left panels in the fig- 
ure show our results without any arbitrary rescaling, while 
in the right panels we apply the empirical rescaling they 
suggest, which effectively removes the trend with mass. 

Park, Kim & Park (2010) have argued that the main 
reason for this discrepancy may be due to the halo def- 
inition itself, and to the halo environment. They used a 
more sophisticated halo finding algorithm (from Kim & Park 
2006) which identifies gravitationally self-bound and tidally 
stable halos. They showed that the dependence of A1/A3 
on local density is stronger for more massive isolated ha- 
los, and argued that tidal interactions between neighboring 
halos make them more spherical on average. In particular, 
Park, Kim & Park (2010) pointed out that the main prob- 
lem relies primarily on the FOF halo finder. This is because 
the FOF identification scheme cannot distinguish between 
objects connected by filaments. For example, if two mas- 
sive and almost spherical neighboring halos are slightly con- 
nected by a filament, the FOF halo finder will consider those 
objects as one big structure, and erroneously assign a rather 
elongated shape to it. Therefore, the FOF scheme will find 
that more massive objects are more elongated, while in re- 
ality this is an artifact of the FOF procedure. Park et al. 
(2010) went one step further; after running their halo finder 
algorithm, they identified subhalos within FOF halos, and 
classified their FOF halos into central (i.e. the most massive 
halos in each group of halos and subhalos), satellite (i.e. 
halos which are not the central ones in each group of ha- 
los and subhalos), and isolated halos (i.e. structures which 
do not overlap with other halos, but which can have close 
neighbor halos and subhalos). After this classification, they 
found that the degree of halo sphericity increases with mass 
(i.e. more massive objects are rounder) for isolated halos, 
up to about 2 x W 14 ^ 1 Mq - a similar trend we find in our 
theoretical study. 

4.2 Conditional axis-ratio distributions 

The left panel of Figure [8] compares the conditional minor- 
to-intermediate axis distribution of A\ /A2 for a given range 
of A1/A3 in our model (histograms), with Jing & Suto's 
simulations (solid lines) at z = 0. In this case, instead, we 
find remarkably good agreement with the numerical mea- 
surements. 

Recall that we had argued that the conditional prolate- 
ness distribution g(p\e, 8, a) in the initial conditions (equa- 
tion [25] and Figure [5}, depends only weakly on 5 /a, so we 
expect the conditional minor-to-intermediate axis ratio dis- 
tribution of the evolved object, p(Ai2| A13, 5, a), to be simi- 
larly universal - independent of mass and time; whether this 
equivalence is due to statistics, rather than physics, is sub- 
ject of ongoing work. The agreement with Jing & Suto's sim- 
ulations suggests that p(Ai2\Ai3,5,a) is indeed preserved 
from the initial conditions. In fact, except for the exponen- 
tial factor in our equation (|25[). our theoretically motivated 
expression for this distribution is exactly the same as that 
of Jing & Suto's empirical fitting formula. 
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Figure 7. Distribution of the minor-to-major axis ratio, A\/A 3 , in our model (dotted) and in Jing & Suto's (2002) numerical simulations 
(solid). Raw data are plotted in the left panels, while in the right panels data are rescaled according to their prescription. 



More recently, Wang et al. (2010) have also found that, 
in their simulations, p(^4i2|-4i3) is independent of mass and 
epoch. They also noted that the final redshift dependence 
of the short-to-intermediate axial ratio is much weaker than 
the other two ratios, indicating that new material tends to 
be accreted along the major axes of halos. Our formula (|25[) 
provides the theoretical/physical explanation for their find- 
ings. 

In this context, we note that Lee, Jing & Suto (2005) 
also presented an analytic expression for the axis ratio 
distribution of triaxial objects, and claimed that it suc- 
cessfully reproduces the conditional intermediate-to-major 
axis ratio distribution. However, in their Figure 5 which 
is supposed to support this claim, it is the minor-to- 
intermediate distribution from the numerical simulations, 
p(Ai/A 2 \Ai/A 3 ), which is compared with their theory 
curve for the intermediate-to-major axis ratio distribution, 
^(^/^alAi/As). While this invalidates their claim (see also 
Appendix [Alfor more discussion on the flaw in their model), 
it does raise the question of why the formula for one distri- 
bution happened to describe another? 



Our model shows why. When Ai <C 1, then A1/A2 — 
1 - <J(ti)(e + p), while A2/A3 ~ 1 - <5(ti)(e - p). How- 
ever, the sign difference in the prolateness does not al- 
ter the Jacobian of the two different transformations: 
(e,p -> A 2 /A 3 ,A 1 /A 3 ) and (e,p -> A1/A2, A1/A3). Hence, 
p{A 3 /A 3 \A 1 /A 3 ) d(A 2 /A 3 ) = P (A 1 /A 2 \A 1 /A 3 )d(A 1 /A 2 ), 
and therefore 

P(A 2 /A :i \ A1/A3) = ■ P{Ai/A 2 \A 1 /A 3 ). (27) 

(A 2 /A 3 ) 

Thus, at fixed A1/A3, the minor-to-intermediate or the 
intermediate-to-major axis ratios distributions are both 
equivalent, and so 'universal'. Although equivalent, clearly 
these distributions are indeed different and neatly distin- 
guishable: the right panel of Figure [8] compares now the con- 
ditional intermediate-to-major axis distribution of A 2 /A 3 
for a given range of ^4i/j43 in our model (histograms), with 
Jing & Suto's simulations (solid lines) at z = 0. As evident 
from the plots, p(A 2 /A 3 \Ai/A 3 ) is more asymmetric than 
p(Ai/^42|^4i/^4s), and the difference is more significant for 
lower values of Ai/A 3 . 



Dark matter halo shapes 11 





0.6 0.8 



Figure 8. Conditional distributions (A1/A2IA1/A3) - left panel - or (A2 / A$\Ai / A3) - right panel - in our model (histograms), and in 
the numerical simulations of Jing &; Suto (2002) (solid lines). Results are shown for two different mass ranges (both at z = 0) and for 
different values of A1/A3, as labeled. 



4.3 Caveats 

In general, a direct comparison between our theoretical 
model and iV-body simulations is both challenging and sub- 
tle. This is primarily because there is still no unanimous 
agreement about how to define a gravitationally bound ob- 
ject; i.e the very definition of what a non-spherical dark 
matter halo is, both in simulations and observationally, is 
a non-trivial problem. Determining which material belongs 
to a halo and what lies beyond is an open problem (Prada 
et al. 2006; Bett et al. 2007; Diemand, Kuhlen & Madau 
2007; Cuesta et al. 2008). Common methods for defining 
virialized halos in simulations include: spherical overdensity, 
tree algorithms based on the branches of the halo merger 
trees, two-step procedures with post-processing, maximum 
circular velocity, etc., (see Klypin et al. 2010, and references 
therein for more discussion). Clearly, the distribution of halo 
shapes measured from simulations depends critically on the 
halo definition. For instance, Bett et al. (2007) found that 
spherical overdensity halos are more spherical than FOF or 
TREE halos, and that FOF halos show a much broader dis- 
tribution of shapes and a strong preference for prolateness. 
This fact complicates the comparison and interpretation of 
theory against any numerical study. 

Since we are most interested in halo shapes, measure- 
ments which do not perform spherical averages are more 
closely related to our models. Jing & Suto (2002) used one 
such method: a friend-of-friend algorithm (Davis et al. 1985; 



Lacey & Cole 1994). This is a percolation scheme that links 
together all particles that are closer than b = 0.2d, where d 
is the mean interparticle separation. This value of b returns 
objects which are approximately 200 times the background 
density: it makes no assumptions about the shape of the re- 
sulting object. However, it may group distinct halos together 
into the same object, confusing the comparison with theory 
(White 2001; Tinker et al. 2008; Lukic et al. 2009). This 
problem is more serious for halos in higher density regions; 
since high density regions have top-heavy mass functions 
(Sheth & Tormen 2002), massive halos identified by this al- 
gorithm may be erroneously considered to be elongated. 

Jing & Suto (2002) used 6 = O.ld, so their FOF clumps 
are smaller and denser than those defined using the more 
conventional 0.2d. As a result, they are less likely to suf- 
fer from spuriously linked halos. However, in effect, their 
procedure identifies only the central parts of the more con- 
ventional halo. 

As a further complication, having settled on a halo def- 
inition, some authors only study shapes of a subset of the 
halo population. E.g., Jing & Suto (2002), only analyzed ha- 
los which looked 'relaxed': at fixed mass, this is almost cer- 
tainly a more homogeneous subset of the entire population. 
Our model does not account for this additional selection. 

Despite such difficulties in the model/simulation com- 
parison, we believe such an exercise is still useful, because 
there are other reasons why we anticipate disagreement be- 
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tween model and simulations. E.g., the ellipsoidal collapse 
itself is a significant oversimplification of the more complex 
nonlinear gravitational evolution - it cannot be expected to 
fully capture the dynamics of iV-body simulations. Despite 
all these caveats, we believe the remarkably good agreement 
between the measured conditional axis ratio distributions 
and our equation (f25[l (Figure [8]) is non-trivial. 



5 DISCUSSION 

We presented a model to describe the shape distribution of 
bound objects based on Rossi (2008), which is a simple ex- 
tension of that introduced by Bond & Myers (f 996) and used 
by Sheth, Mo & Tormen (2001) to estimate the abundance 
by mass of these objects. Our analytic prescription is made 
of two independent parts: one is a scheme for how an initially 
spherical patch evolves and virializes, for which we adopted 
the ellipsoidal collapse dynamics (Section [2l Figure [l}; the 
other is the correct assignment of initial shapes to halos of 
different masses, achieved through the excursion set formal- 
ism (Section [3] Figure |4|. Along the way, we discussed an 
analytic approximation to the evolution, considerably more 
rigorous than the Zeldovich approximation fSection !2.2l Fig- 
ure [2j. And we introduced a useful planar representation of 
halo axis ratios, which gives a simple intuitive framework for 
discussing the dynamics of collapsing regions (Section 12.31 
Figure [3j). 

We showed that the model is able to provide a rea- 
sonable description of the minor-to-major axis ratio dis- 
tribution seen by Jing & Suto (2002) in their numerical 
simulations only within a limited mass range (i.e. around 
M ~ 10 13 /i _1 Af Q ~ M„). Outside this range, Jing & Suto 
found that low-mass halos are preferentially more spherical 
than high-mass halos, whereas in our model, the high-mass 
halos are more spherical (Figures [15] and 0). We argued that 
some the disagreements may originate from differences in 
halo types and mass range, and may be alleviated if we con- 
sider only isolated halos. This is because Park, Kim & Park 
(2010) find that, on average, high-mass isolated halos are in- 
deed more spherical than low-mass ones, and the difference 
is larger in high density /shear regions (see their Section 3.4 
and their Figures 9, 10 and 11, for details on how to select 
isolated clusters from simulations). Note also that isolated 
clusters can be readily identified from an observational cat- 
alog, as done for example in Park & Hwang (2009) for a 
sample of Abell galaxy clusters (their Section 2.2 gives a de- 
tailed explanation on how the selection is performed). The 
finding of Park, Kim & Park (2010) has an interesting con- 
nection to recent observational work on mass scales where 
the hypothesis that a halo hosts only one galaxy may be 
appropriate. For early-type galaxies in the SDSS, although 
the observed projected axis ratio b/a increases with increas- 
ing stellar mass or luminosity, it decreases at the highest 
masses (Bernardi et al. 2008, 2010). This reversal is thought 
to arise because the most massive galaxies tend to be in re- 
gions (e.g., at cluster centers) where recent radial mergers 
may have made them prolate. Fasano et al. (2010) come to a 
similar conclusion based on an analysis of the intrinsic shape 
distribution of BCGs from the WINGS survey (Fasano et al. 
2006; Varela et al. 2009). They propose that the prolateness 
of the BCGs (in particular the cDs) could reflect the shape 



of the associated dark matter halos (also see Smith et al. 
2010). Hence, while it is tempting to conclude that the dis- 
crepancy between mass dependent trends in simulations and 
theory will be alleviated if we only consider isolated halos, 
we note that Jing & Suto (2002) eliminated halos in their 
simulations which were not relaxed. If these had suffered re- 
cent mergers, then they may already have removed the most 
prolate halos from their sample. Further investigation of this 
point is the subject of ongoing work. 

Despite this disagreement about the distribution of 
minor-to- major axis ratio, we found very good agreement be- 
tween our model and the conditional minor-to-intermediate 
axis ratio distribution measured from simulations (Figure 
[5|. In particular, we showed that our model provides phys- 
ical motivation for the empirical fitting formula obtained 
in Jing & Suto (2002) from numerical studies (our equa- 
tion [25}. In our model, this distribution is closely related to 
the conditional prolateness distribution in the initial condi- 
tions (Section 13.21 Figure [5}. Unfortunately, results in Lam 
et al. (2009) suggest that this distribution is unlikely to be 
able to discriminate between Gaussian initial conditions and 
ones where the primordial non-Gaussianity is of the 'local' 
type (i.e. Rossi, Chingangbam & Park 2010). 

We also discussed the agreement /disagreement between 
theory and numerical predictions, in the context of known 
problems and limitations in the modelling to difficulties in 
making shape measurements from simulations. 

We are currently studying the following improvements 
and extensions to our model: 

(i) Our model makes the extremely simple assumption 
that axis lengths freeze-out once they have shrunk to a suf- 
ficiently small size. Subsequent violent relaxation effects as- 
sociated with particles which collapse along other axes may 
change this size - something our current model does not 
consider. For instance, massive halos are expected to have 
assembled their mass more recently than low mass halos. 
Therefore, relaxation effects may have had more time to 
change the axis ratios of low mass halos than of higher mass 
halos. If the net effect of this relaxation is to make the halos 
more spherical than they would otherwise have been, then 
accounting for it may help resolve the discrepancy between 
our model and simulations. Note that halos having smaller 
amounts of substructures tend to be closer to virial equilib- 
rium (e.g. Giocoli et al. 2010), but there has been no study of 
whether or not haloes with abnormally small sub-structure 
(for their mass) are rounder. 

(ii) Testing different collapse criteria, and other prescrip- 
tions for the external tidal field (e.g., Angrick & Bartelmann 
2010). 

(iii) Accounting for correlations between the properties of 
halos and their environment (assembly bias), as seen in N- 
body simulations (Sheth & Tormen 2004; Ragone-Figueroa 
& Plionis 2007) . Environmental effects can have a significant 
impact on the properties of virialized halos because, in the 
ellipsoidal model, the critical density threshold <5 CC depends 
strongly on the initial values of e and p, and these distri- 
butions depend on environment (Keselman & Nusser 2007; 
Wang, Mo & Jing 2007; Desjacques 2008). In our current 
implementation, only the initial density of the surrounding 
environment matters. 

(iv) Accounting for correlations between formation times 
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and halo shapes, as seen in TV-body simulations (Ragone- 
Figueroa et al. 2010). 

(v) Modifying the algorithm of Sheth & Tormen (2002) 
to account for correlated steps, following Maggiore & Riotto 
(2010) and De Simone et al. (2010), when deriving our initial 
conditions. 

(vi) Including the effects of baryonic physics. The con- 
densation of baryons to the centers of dark matter halos 
is known to change their shape (Dubinski 1994; Holley- 
Bockelmann et al. 2002; Kazantzidis et al. 2004; Debattista 
et al. 2008). 

(vii) Estimating correlations between halo shapes. This is 
possible because, in our model, the final shape is determined 
by the initial conditions, and correlations in the initial con- 
ditions can be quantified. Thus, in our model, correlations 
at the present time are obtained by taking the appropri- 
ate average over the initial correlations - this average being 
determined by the mapping from (e,p) — > (A2/A3, A1/A3). 
Most published estimates of the correlation between shapes 
do not account for this mapping. 

(viii) Describing the shapes of subhalos as well as voids 
(e.g., Shandarin et al. 2006; D'Aloisio & Furlanetto 2007; 
Platen et al. 2008). 

We conclude with the observation that while numeri- 
cal studies are indispensable for quantifying the shapes and 
structures of dark matter halos, we hope that our simple 
procedure may provide a complementary theoretical frame- 
work for understanding halo shapes, and how these shapes 
are correlated with larger-scale structures. 
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APPENDIX A: COMPARISON WITH THE 
MODEL OF LEE, JING & SUTO (2005) 

The main text describes a number of inconsistencies in the 
work of Lee et al. (2005). Besides the mismatch between 
p(A 1 /A 2 \A 1 /A 3 ) and p(A2/A 3 \A 1 /A 3 ) (see Sectiong^, we 
also stated in the main text that their approach does not 
reduce to the Zeldovich approximation - despite the fact 
that they believe it does. 

To see why, note that their equation (2) gives the usual 
expression for the nonlinear density in the Zeldovich ap- 
proximation - the same one we use. Namely, if we use Ri 
to denote the initial length of axis i, and n its length at 
a later time, then their equation (2) comes from setting 
ri/Ri = 1 — D(t)Xi, where Ai is the eigenvalue of the ini- 
tial shear tensor. Note that there is no square root factor 
in the relation between Ai and the Zeldovich-evolved axis 
length ri/Ri. 



Lee et al. (2005) then define (a, b, c) to be the eigenval- 
ues of the inertia tensor. For a homogeneous ellipsoid with 
axis ratios (A,B,C), the inertia tensor eigenvalues (a,b,c) 
are related to (A 2 , B 2 , C 2 ), whereas their equation (7) sets 
a equal to the square root of A. I.e., they have put the 
square root on the wrong side of the equation. Therefore, 
their expressions will not reduce correctly to the Zeldovich 
approximation at early times. 

Note that their equation (8) does in fact have the cor- 
rect relation between the semi-axis lengths and eigenvalues 
(because it has been taken correctly from Equation 7.2 of 
Bardeen et al. 1986). In this notation, (i = a 2 where the 
£i are the eigenvalues of the inertia tensor, but Lee et al. 
(2005) incorrectly call (a, b, c) the eigenvalues of inertial ten- 
sor, when they are in fact the semi-axis lengths (see text 
above Equation 7.2 of Bardeen et al. 1986). However, by 
equation (15), Lee et al. are referring to b/c as ratios of axis 
lengths, rather than as ratios of eigenvalues. Their sloppiness 
about the important difference between these two gives rise 
to the following problem. If (a, b, c) are axis lengths, then the 
square root factor in their equation (7) means they will not 
recover the Zeldovich approximation, which does not have a 
square root factor. If they are eigenvalues of the inertia ten- 
sor, then they have put the square root in the wrong place, 
so again, they will not recover the Zeldovich approximation. 

Unfortunately, Jing & Suto (2002) are also ambiguous 
about just what it is they measured. Their equation (5) 
shows (a, 6, c) as being axis ratios, whereas their text de- 
scribes them as being obtained from diagonalizing the iner- 
tia tensor, hence eigenvalues (unfortunately, they don't say 
explicitly which is the square-root of the other). 

In a subsequent study on cosmic voids, Lavaux & Wan- 
delt (2010) used some erroneous definitions from Lee et al. 
(2005) but still found excellent agreement between their an- 
alytic framework and measurements from numerical simula- 
tions. This is because they where "consistently wrong", in 
that their derived quantities measured from the simulations 
suffer from the same flaw present in their theory. Hence, 
although their ellipticity definitions are incorrect, the com- 
parison theory/simulations is still fair. We also note that the 
philosophy of our approach in modelling halo shapes is simi- 
lar to theirs, although the details are different (see their Ap- 
pendix B) and they matter (i.e. in particular the use of the 
ellipsoidal collapse barrier as in Equation 1221 which allows 
one to distinguish between peaks and random positions). We 
will expand more on this important issue in a forthcoming 
publication. 



